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Abstract 

We have performed two-dimensional lattice-gas-automaton simulations of 
granular flow between two parallel planes. We find that the velocity profiles 
have non-parabolic distributions while simultaneously the density profiles are 
non-uniform. Under non-slip boundary conditions, deviation of velocity pro- 
files from the parabolic form of newtonian fluids is found to be characterized 
solely by ratio of maximal velocity at the center to the average velocity, though 
the ratio depends on the model parameters in a complex manner. We also 
find that the maximal velocity {umax) at the center is a linear function of the 
driving force (g) as Umax = ctg — S with non-zero 5 in contrast with newtonian 
fluids. Regarding density profiles, we observe that densities near the bound- 
aries are higher than those in the center. The width of higher densities (above 
the average density) relative to the channel width is a decreasing function of a 
variable which scales with the driving force (g), energy dissipation parameter 
(e) and the width of the system (L) as g^L'^ /e with exponents = 1.4 it 0.1 
and z^ = 0.5ib0.1. A phenomenological theory based on a scaling argument 
is presented to interpret these findings. 
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I. Introduction 



Granular materials exihibt many interesting phenomena. The segregation of par- 
ticles [|l], D, H, ^, heap formation and convection cells under vibration [0, 0, H, 
and anomalous sound propagation ||^ are just a few examples. Such phenomena 
occur because the dynamical responses of granular media are quite complex and 



different from those of usual solids, fluids and gases |10, 11, 12, 13, 14]. Consid- 



ering the complexity in the dynamics, one is tempted to study first the granular 
materials in simple geometries and proceed then to more complicated situations. 
However, even in the simplest geometries such as hoppers and tubes, their flow un- 
der uniform external driving force (such as gravity) still shows complex dynamics 
T3| , |TB|, |T^, [T^, |T^, ^ For tubes, one can observe density waves in granualr 



flow when the width of the tubes is narrow enough. Here the friction of the wall, as 
well as the dissipation among the grains themselve, plays an important role in the 
occurence of density waves. 



The purpose of this paper is to study granular flow in channels (different from 
tubes with wider width) and compare it with flow of Newtonian fluids. Compar- 
ison between granular flow and ordinary fluids is useful because it links with the 
knowledge that we have gained before. Besides the obvious contrast between the 
discretness of granular materials and continuum of fluids, a key difference is that 
collisions of particles in granular flow are inelastic (energy-dissipative). In fluids, 
channel flow is a typical and fundamental topics in rheology (Poiseuille flow). Un- 
der the non-slip boundary conditions valid for viscous fluids, Newtonian fluids have 
a parabolic velocity profile across the channel. To claim that granular flow is non- 
Newtonian flow, a non-parabolic velocity profile is a clear and fundamental evidence. 
In this paper, we will report the non-parabolic shape of velocity profiles and non- 
uniformly distributed density profiles in granular channel flow. 



Study on granular materials has a long history in engineering. Much engineering 
literature is devoted to understanding how to deal with these materials. However, 
it is only after the pioneering work of Bagnold that theoretical work has begun 
to progress. Instead of focusing on the detailed dynamics of individual particles, 



granular hydrodynamics [p3| , ^ treats the granular media like "fluids" . One defines 
a set of macroscopic quantities like particle density, velocity, and granular tempera- 
ture. By assuming local equilibrium |^ one can write down the equations based on 
mass conservation, momentum conservation and energy balance. These equations 
are similar to the Navier-Stokes equations for fiuids. 



Using the idea of hydrodynamics. Savage p6| has studied the channel fiow of 



cohesionless granular materials and obtained theoretically the velocity profiles. He 
obtained the velocity profiles under non-slip boundary conditions for different pa- 
rameter R (which depends on the system parameters and on the density at the 
channel center). In the limiting case of i? = (which corresponds to the situation 
where density are uniformly distributed), the non-dimensional velocity profile has 
the form of f = 1 — x^/"^ where v is the velocity normalized to 1 with respect to the 
maximal velocity at the center, x is the distance from the channel center divided 
by the half-width of the channel. For larger i?, the velocity profiles become more 



blunted (but without simple analytical form). Savage p6[ performed also experi- 
ments to measure the velocity profile and found that his experiments corresponded 
to the case of R = 0.3. Since he did not intend to compare his results with those for 
fiuids, his theory did not recover the parabolic velocity profile in any limiting case 
as he included from the begining that the viscosity is proportional to the shear rate 



1 2^ . However, we find from his experimental data that the velocity is very close to 



have a parabolic shape. 

In the following we will report our numerical results on the velocity profiles. We 
find that when the dissipation is very weak the velocity profiles are close to the 
parabolic shape while in the other cases they are more blunted than the parabolic 
form. As mentioned above, we realize that the key difference between fluids and 
granular flow is that the latter is energy dissipative. Therefore, when the parameter 
controlling the energy dissipation is switched off, the theory should recover that of 
fluids. Furthermore, we observe that density near the boundaries are higher than 
those at the channel center. Our numerical flndings are consistent with the earlier 
experimental observations by Soo et al |^ : "the nature of the concentration, mass 



flow, and velocity distributions of solid particles is such that the concentration in- 



creases toward the wall of the pipe, mass flow decreases toward the wall, and velocity 
is less than or equal to that of the stream at the core" . 



In this paper we study the channel flow by computer simulations. Since a general 
theory for granular media is not yet available, people have used various compu- 
tational methods to get better understanding about the complicated rheological 
behavior of granular media. Among the different methods are molecular dynamics 
(MD) 0, |18|, H, Monte Carlo simulations |, |3|, the diffusing void model 13, 
event driven algorithms and cellular automaton |3^. So far the most widely 
used method is MD which simulates the granular materials on a "microscopic" 



level (the grain's level). MD has been recognized to be very successful in simulating 
granular materials |^ . MD needs, however, much computer time to give reasonable 
results. To calculate the velocity and density profiles one needs long time average 
in order to get reasonable statistics. 



The same situation was also faced in classical fluid mechanics some years ago 
when Frisch, Hasslacher and Pomeau ||3^ proposed lattice-gas automata (LGA) 
as a novel alternative to the direct solution of the equation of motion. As a sort 
of primitive molecular-dynamics system LGA offers the advantage of guaranteed 
numerical stability coupled with extreme computational simplicity. The basic idea 
behind LGA is that a properly defined cellular automaton with appropriate conser- 
vation laws should lead to the Navier-Stokes equations. A detailed comparison of 
simulational results from LGA and the well-established theory of Newtonian fluids 
in Poiseuille flow can be found in Refs. The main check to the LGA for the 

Newtonian dynamics was the parabolic velocity profile in channel flows | 



39, 40 



In this paper we employ the LGA model of Ref. PU| to study the velocity and 
density profiles in granular channel flow. The LGA of Ref. |^0|] is an extension of 
the LGA of usual fluids by including energy dissipation among the particles into the 
model. As emphasized above, energy dissipation is a major mechanism by which 
granular flow differs from ordinary fluids. The LGA models for granular materials 



were sussessful in simulating granualr flow [20, 41, H2, 43 . 



This paper is organized as follows. We describe the model in some detail in Sect. 
II. The results obtained by simulations are presented in Sect. III. A phenomenologi- 
cal argument is made in Sect. IV to understand the numerical results. A discussion 
is contained in Sec. V. 

II. Simulational model 

We consider an LGA at integer time steps i = 0, 1, 2, • • • with N particles located 
at the sites of a two-dimensional triangular lattice. The arrangement of the triangle 
lattice is depicted in Fig. 1 (a). At each site there are seven Boolean states which 
refer to the velocities, Ck{k — 0, 1, 2, • • • 6). Here Ck{k — 1, 2, • • • 6) are the nearest 
neighboring (NN) lattice vectors and cq — refers to the rest (zero velocity) state. 
Each state can be either empty or occupied by a single particle. Therefore, the 
number of particles per site has a maximal value of 7 and a minimal value of 0. In 
this paper we use the number of particles per site as density which can therefore 
be greater than 1.0. The time evolution of the LGA consists of a coUision step 
and a propagation step. In the collision step particles change their velocities due to 
collisions and in the subsequent propagation step particles move in the directions of 
their velocities to the NN sites where they collide again. 

The system is updated in parallel. Only the specified collisions shown in Fig. 1 
can deviate the trajectories of particles. All collisions conserve mass and momentum. 

For two- and three-body collisions, we have the probabilistic rules shown in Fig. 
1 (b). The probablity that a configuration may take place is shown next to the 
configuration. If the parameter e is nonzero, it means that energy can be dissipated 
in the collision. 

Collisions with rest particles may produce more than one rest particles on that 
site. This is allowed temporarily as in the Fig. 1 (c). However, immediately after 
the coUision step, the extra rest particles randomly hop to NN sites until they find 
a site with no rest particle and there they stop. 



We use no-slip boundary conditions at the channel walls which are parallel to 
the y-axis and periodic boundary conditions along the channel. No-slip conditions 
are employed by allowing that any particle colliding with the wall along any of the 
three possible directions bounces back into the incoming direction. 

We incorporate the driving force, namely, gravity, in the way that Kadanoff et al 
Il38| have used. The direction of gravity is along the y-axis, i.e., the downward in 
Fig. la. After each time step we randomly select a lattice site and, if possible, apply 
one of the forcing rules: (i) a rest particle goes into motion with equal probability 
along one of the two lattice directions which form an angle of 30° with the direction 
of gravity; (ii) a moving particle change its velocity by a unit vector along the direc- 
tion of gravity if the resulting vector is possible on the triangle lattice used. Each 
successful application of a forcing rule adds one unit of momentum to the system. 
The forcing process is repeated until the desired amount of momentum (we label it 
as G) has been transfered to the system; fractional amounts of momentum to be 
added to the system are accumulated across time steps until they sum to an amount 
greater than 1, at which time one additional unit of momentum is added to the gas. 
The actual forcing scheme is slightly more complicated since it must compensate 
for inhomogeneity in the momentum and number densities due to the macrosopic 
flow (see also Ref. ||38[). The forcing algorithm randomly selects a lattice row and 
column and then searches along that column until it finds a site where a forcing 
rule may be successfully applied. This guarantees that forcing operations will be 
uniformly distributed across the width of the channel, despite variations in the mass 
and momentum densities. It is noted that the magnitude of the gravity cannot be 
extremely large in our model system due to the following fact: on one hand the 
maximal velocity at the channel center increases with the increasing gravity, on the 
other hand the LGA only guarantees a velocity less than unit. 

III. Simulational Results 

We evolve the system according to the collision rules defined above. The initial 
configuration of the system is set to be random in the sense that every state (except 
the rest state) of each site is randomly occupied according to a preassigned average 



density p. We discard the configurations in the first period of several thousand time 
steps until we are sure that the system is in the steady state where energy input by 
gravity is averagely equal to the energy dissipation in the system. Fig. 2 shows a 
typical curve of the total kinetic energy relaxation starting from a random configu- 
ration at time step t = 0. 



After the kinetic energy curve becomes flat, we then make the necessary averge 
over space and time as follows. Let us label the particle number of state k{k — 
0, 1, 2, ■ ■ ■ , 6) on the ith lattice site with coordinates (xj, i/i) at time step t as nk{xi, i/i, t). 
As depicted in Fig. 1(a), the x-coordinates of lattice sites take values of inte- 
gers and half-integers of 0, 1/2, 1, • • • , L — 1, while the y-coordinates have values of 
0, 1/2, 1, • • • , M-1 with M even. The total number of lattice sites is N ^ LM-M/2. 
In this paper, we always keep M — 2L and fix the average density p = 1.5 unless 
otherwise specified. 

We are interested in the density profile p{x) and velocity profile in y-direction 
u{x) — Vy{x) with X integer 

^ to+T N 6 
^-^ -'^^ t=to+l i=l k=0 

to+T N 6 

^ OTM E Y.J2(^k^MXi,yi,t)S{[Xi]-x) (2) 
^-^ t=to+l i=l k=0 

where the function [z] takes the integer value of its variable z and function 6{z) equals 
to one when z—0 and zero otherwise. Ck — {Ck \ cj^ ^) = {cos{7r{k — l)/3), sin{7r{k — 
l)/3)) for (k — 1, 2, • • • , 6) are the unit vectors of the lattice bond orientation (and 
Co — 6). We actually make average over two column whose x-coordinates have the 
same integer part, they are columns with Xi (integer) and Xi + 1/2 (half-integer), 
therefore, the normalization factor is where T is the time period over which we 
make average starting from time step Iq. T is usually several tens of thousand time 
steps. 



As a first step, we check that when the dissipation parameter e is set to zero 
we recover the parabolic velocity profile of Newtonian fluids. This is shown in Fig. 
3(a) for system size L — 64 with average density p — 1.8. From Fig. 3(a) we can 



estimste the magnitude of the kinematic viscosity equal to 0.35 according to Ref. 
||38|| . Meanwhile, the density profile in Fig. 3(b) is just a flat curve, showing that 
density is uniformly distributed. 

Fig. 4(a) and Fig. 4(b) show the velocity and density profiles when energy dis- 
sipation with e = 0.01 is present in the system. We see clearly that the parabolic 
velocity profile is no longer kept in this case, but more blunted. In fact, a parabolic 
curve demanding the maximum value at the center equal to the numerial data and 
zero at the two boundaries deviates substantially from the numerical data. It is also 
remarked that the density is not uniformly distributed, in contrast with those of 
Newtonian fluids. 

Regarding the velocity profile, we can generally write down the following equation: 

u{x) = u^ax{e,g,L)u{^,eL''\gL''^) (3) 

where Umax is the maximum velocity. Here gravity is taken to be system size inde- 
pendent quantity g = G/N where G is the total amount of momentum added to the 
system after each time step and is the total number of lattice sites. 

A simple measure of the flatness of velocity profiles is the ratio 7 = Umax I u where 
u is the mean value of u{x). It is equal to 1.5 for Poiseuille flow of Newtonian fluids 
and 1.0 for a perfectly flat profile. Note that 7 can be written as 

= [\{x,eL''\gL''^)dx (4) 
Jo 

Fig. 5 displays the ratio 7 versus gravity for two different dissipation parameters 
e = 0.01 and e = 0.02 with different channel widths. From Fig. 5 we see 7 depends 
on the channel width, degress of dissipation and gravity. We find that 7 can be solely 
used to characterize the shape function u{x^eL"^ , gL""^). Taking two set of different 
parametres which leads to almost identical 7, e.g., {ei = 0.01,(71 = 8.0, Li = 64} 
and {t2 = 0.02,^2 = 10.0,^2 = 48} resulting 71 = 1.27594 and 72 = 1.28298, 
we show in Fig. 6 that their shape functions are almost identical. Note that both 
set of parameters {e,g,L} share no common value but their resulting 7 are almost 
identical. For comparison, the velocity profile for a different value of 7 = 1.35906 



is also plotted in Fig. 6. As expressed in Eq. @), 7"^ is the 0th moment of the 
distribution function u{x). Fig. 6 tells us that the distribution function u has such 
a property that as long as the 0th moment is determined higher order moments are 
determined. This may suggest that the parameters {e, g, L} combine to give a scalar 
on which u{x) depends solely. 

In order to fit the velocity data, we have tried the polynominal fitting routines. 
We find that polynominal up to an order of 4 is not sufficient; most of the data can 
be fitted by a polynominal up to an order of 6: 

u{x) = Vq — V2{x — cY — Vi{x — cY — Vq{x — c)^ (5) 

where x is the normalized coordinate in the interval [0:1] and c is the coordinate of 
the channel center. The three velocity profiles are fitted using Eq. (|^) in Fig. 6. 
In the next section basing on a phenomenological theory we will propose a formula 
for velocity profiles where the dependences of these fitting parameters on the model 
parameters {e,g,L} are clearly expressed and we will check our theory by the sim- 
ulational data. 

For Umaxi^, 9, L), we find that it is a hnear function of g as Umaxi^, 9, L) = ag — S 
with non-zero 6. Fig. 7 displays four curves of Umax for different channel widths but 
with fixed e. The slope a depends on L but not as which is valid for Newtonian 
fiuids. If we increase e, the value of a decreases. Non-zero 6 means that the driving 
force must exceed a non-zero threshold {6 /a) in order to drive the granular media 
into motion. From the statistics we are not sure whether 6 is independent of L 
but one can see from Fig. 7 that the depenence is very weak, if any. Futhermore, 
although not shown in Figure, we find that 6 have a tendency to approach zero as 
e approaches zero, which recovers the case of Newtonian fiuids. 

Regarding the density profile, we observe that densities near the walls are higher 
than those at the channel center, in contrast to the uniform density distribution of 
Newtonian fiuids. A simple way to characterize the spatial variation is to define 
a width over which densities are higher than the average density. We plot such a 
width (A) divided by channel width in Fig. 8 (a) against gravity g for different 



system sizes and different dissipation parameters. Generally, A/L decreases with 
increasing g and increasing L but increases with increaing e. The scattered data in 
Fig. 8(a) can be scaled to one curve as shown in Fig. 8 (b) where we plot A/L 
versus g^L'^ /e. We find that the exponents are /x = 1.4 ±0.1 and v = 0.5 ± 0.1. 



We have checked that the above observations do not change their properties qual- 
itatively as we change the average density p. The scaling exponents fi and z/ are 
independent of the average density. 



IV. Phenomenological Theory 



In this section we present a phenomenological theory to interpret our numerical 
findings. We note that there are two important system parameters, namely, g and 
e: g is the velocity change per particle per time step and e is energy dissipation per 
particle per time step. Therefore, using dimension analyse, we have 

dV dE 

where V and E are characteristic velocity and kientic energy of a particle. Here- 
inafter we put the particle's mass to be unity. Since V"^ ~ -E, gV/t must be dimen- 
sionless. 



We may construct a quantity which has a dimension of velocity by using the 
channel width L as 

V=4^ (7) 
The above argument implies that P defined as 

^ gV g^'^L^/^ 

P = — = (8) 

must be a fundamental dimensionless quantity. Eq. (|^) is consistent with the sim- 
ulations illustrated in Fig. 8b where A/L is a function of a variable which scales as 
g^L'^ /e with exponents p = 1.4 ± 0.1 and u = 0.5 ± 0.1. 



Note that when e = 0, the system has another dimensionless quantity which can 



be derived from the newtonian hydrodynamics equation 

g + = (9) 

where z/ = /i/p is the kinematic viscosity. The two terms in Eq. have the same 
dimension, therefore, 

Q = -T7 = = (^0) 



is another dimensionless quantity. Note in Eq. (|T0|) we have used Eq. (|^). Q is in 
fact the square root of the Reynolds number. 

Now, we consider the maximum velocity u^ax at the channel center. Simulations 
show that Umax is a linear function of g 

Umax = ag - S. (11) 

We determine the forms of a and 6. First of all, they do not depend on g by 
definition. Using Eq. (^, a must take the form 

a = ^«(P,Q) (12) 

The form of a can be determined by the following conditions: (i) a (instead of 
a) is independent of g; (ii) when e (P ^ oo), Umax = gL"^/^- The simplest 
interpolation formula for a is given by 

where ci,C2 are dimensionless positive numerical constants. In fact we have 

1 



c^g = ygL 



1/3 



g3/2ll/2j ^gl/2l3/2 

gVL 



(14) 



e \V3 u 



Therefore, 



L . . 

Cl (eL) ' + C2U/L 



When e = 0, we recover the well-known result: 



Umax ~ — (16) 



When e ^ oo, Eq. (|T^) shows 



a 



cie 



1/3 



(17) 



Our numerical results also show that a decreases with increasing e, in the trend 
of Eq. (0). The L^/^ is even reached in the simulations with e = 0.02. This is 



illustrated in Fig. 9. From Eq. (15) one may see that as long as e is non-zero, the 
term with ci dominates the denominator for large L, leading to the validity of Eq. 



Now we determine 5. It should vanish at e = and be independent of g for any 
non-zero value of e. Using again the fact of Eq. (|^) 6 should be given by 

S = c,^p-'/'f{P/Q') = C3(eL)V3/(p/Q3) (18) 

where C3 is a positive constant. Note that P/Q^ = u^/i^eL^) is independent of 
g. It seems that the form of f{x) can not be specified phenomenologically. If we 
impose the fact that 6 is insensitive to L as observed in simulations, f{x) must be 
f{x) ~ x^^^'^ so that we have 

6 = C3(eL)i/Vi/4e-i/i2L-i/3 = c^ieuy/^ (19) 

Now we turn to the shape function 'u(x) of velocity profile. Since from the sim- 
ulations we know the velocity profile can be well fit using polynominal up to 6th 
order, we propose phenomenologically 

u{x) = 1 - 4u2(x - 1/2)2 _ iQy^^i^x _ 1/2)^ - 6Auq{x - l/2f (20) 
where the coefficients U2, M4, Uq are functions of e, g and L. 

As expressed in Eq. (p, the 0th moment of 'u(x) is characterized by 7. We find 
numerically that the scattered data of 7 in Fig. 5 can be collapsed to one curve by 
using a combination of the dimensionless quantities P and Q: / P = g^^'^L^^^'^e/ . 
Fig. 10 illustrates this point by plotting 7 ag ainst F = IQ-^g^/^L^^/^e for the data 
of Fig. 5. We propose that the coefficients U2, U/^.uq in Eq. (PP]) are functions of F 
and we determine their forms as follows. 



Setting X = in Eq. ( ]20|) and requiring ^(O) = u{l) = 0, we know that U2, U4, Uq 
are related by 

U2 + Ui + Uq = 1 (21) 

So we need only to determine U2 and U4 which must satisfy the following conditions: 

(i) when e = 0, M2 = 1, M4 = which recovers the parabolic velocity profile for fluids; 

(ii) when e —* 00, the velocity profile becomes completely fiat and therefore U2 = 0, 
M4 = O.We propose U2 and M4 as functions of F: 

= rhr '^^> 



where d, e, h are constants. Substituting the above expressions into Eq. ( pUD and 
integrating that equation we obtain 

_^ _ 90 - 2OM2 - Qua 



7 



105 

(24) 



90 



20 6er 



105 

We ckeck the above expressions by fitting the numerical data of 7 using Eq. (0). 
This is shown in Fig. 10. One sees that the fitting is very well except for 7 close to 
1.0 (corresponding to rather fiat velocity profiles where {x — 1/2)^ might be neces- 
sary in Eq. (^0]).). 



Finally we make a remark that velocity profile such that u{x) = l — Ui^{x — l/2)'^ — 
Au2{x — 1/2)^ with < < 2 is not in accord with 7 shown in Fig. 5. 



V. Discussions 



In this paper we have numerically studied the granular fiow in channels using a 
two-dimensional lattice gas automaton model. We have explored the dependences 
of velocity profile on dissipation, gravity and channel width and compared them 
with those of Newtonian fiuids. Parabolic velocity profile is no longer valid as long 
as there is energy dissipation in the system. However, the deviation from parabolic 
profile depends on the degree of energy dissipation in a smooth way (there is no sharp 



transition near e = 0). This may explain the fact that the experimental data of Sav- 
age |2B] is close to a parabolic velocity profile, though he did not notice this point 
in his paper. However, the velocity profile can also be very flat depending on the 
dissipation and the channel width. The observation of Soo confirms this point 
experimentally. Our results also show that the velocity profiles are more blunted 
than the parabolic form. We do not find any velocity profile that is sharper than 
parabolic shape (i.e., our 7 are all less than 1.5), in contrast with Savage's theory 
p6| which predicted that profiles sharper than parabolic ones are possible to exist. 
The maximal velocity at the center is a linear function of gravity as Umax = ag — 5 
but with non-zero 5. Here a behaves no longer as as for Newtonian fluids. Den- 
sity profiles are no longer uniform but are higher near the boundaries than at the 
center. From the simulations we note that this nonuniformity is caused by the rest 
particles while the moving particles are distributed uniformly. This may have a link 
to the clustering of dissipative systems |Q. The boundaries which are represented 
by rough walls seem to serve as clustering seeds due to the non-slip boundary condi- 



tions we used. An earlier experiment by Soo ^Iso showed that density increases 
toward the wall, as we find here. It is interesting to note that the recent experiment 



of Pouliquen and Gutfraind showed that density near the wall is lower than 
that at the center. We attribute this contradiction with our numerical results to 
the different boundary conditions. In their experiments, considerable slip velocities 
were observed at the walls. Both increasing and decreasing density profiles toward 
the walls in different regimes were obtained theoretically in Ref. [^] and both kinds 
of experimental evidences were cited there. 



In our model we find that the width of the region of higher densities (than the 
average density) relative to the channel width is a decreasing function of a variable 
which scales with the driving force ((7), energy dissipation parameter (e) and the 
width of the system (L) as g^U je with exponents yU = 1.4 ± 0.1 and v = 0.5 ± 0.1. 
Such a scaling was not found in any previous investigations. We have shown that 
is a fundametal dimensionless parameter in the present system. 
Combined this with Q = g^^'^L^^'^ /u, all of the results obtained in simulations have 
been successfully understood. 



In this paper we have used non-shp boundary conditions at the channel walls. 
Unlike in viscous fluids, non-slip boundary conditions are not necessarily satisfied 
in every case of granular flow In fact, the general boundary conditions for gran- 
ular flow is still lacking. However, as noted by Savage [^, the non-slip condition 
may be realized by providing wall roughness of the same order as that of the parti- 
cle surface constituting the granular material, as he achieved in his experiments PB 
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Figure Captions 

Figure 1: (a) Sketch of the two-dimensional triangle lattice. Gravity is along 
y-direction. One complete hexagon is drawn; (b) Probabilistic collision rules for 
two- and three-body collisions. Thin arrows represent particles and small circles 
stand for rest particles. The number next to a configuration is the probability that 
the configuration takes place; (c) Collision rules for moving particles with a rest 
particle. Immediately after the collision, more than one rest particle on a site will 
hop to the nearest neighbouring sites randomly until they find a suitable site with 
no rest particle already there. 

Figure 2: Total kinetic energy E versus time step by running the LGA starting 
from a random initial configuration. The kinetic energy of a moving particle is taken 
to be a unit. The plotted curve was obtained with e = 0.01. The total momentum 
added to the system each time is G = 7.0, channel width L — 64 and average density 
p= 1.5. 

Figure 3: (a) Velocity profile for the LGA with e set to zero. Here the total 
momentum added to the system each time is G = 1.0, channel width L = 64 and 
average density p — 1.8. The curve is the best (least-square) fit to the numerical 
data using a polynomial fitting rountine up to second order; (b) Density profile cor- 
responding to (a). 

Figure 4: (a) Velocity profile for the LGA with e = 0.01. Here the total momen- 
tum added to the system each time is G — 7.0, channel width L — 64 and average 
density p — 1.5. The upper curve is the best (least-square) fit to the numerical 
data using a polynomial fitting rountine up to second order and the lower curve is 
a parabolic curve by demanding the maxmial value at the center equal to the nu- 
merial data and zero at the two boundaries; (b) Density profile corresponding to (a). 

Figure 5: The ratio 7 = Umax/u versus gravity g for different channel widths 
L — 32,48,64, 128 (from top to bottom) with e = 0.01 (data points connected by 
solid fines) and e = 0.02 (data points connected by dashed lines). 



Figure 6: The normalized shape function u{x) as a function of x in the inteveral 
[0:1]. Two curves corresponding to {e = 0.02, L = 48, G = 10.0(^ = 0.00434028)} 
{0) and {e = 0.01, L = 64, G = 8.0(^ = 0.00195312)} (+) are almost identi- 
cal, though the two set parameters share no common value. Another curve (□) 
corresponding to {e = 0.02, L = 48, G = 5.0{g = 0.00217014)} is shown for com- 
parison. This curve shares two parameters e and L with the curve {(}) but their 
difference is obvious. Lines are best fit to the numerical data using a polynomial 
fitting rountine up to 6th order as expressed in Eq. (H). The fitting parameters 
are {vq = 0.999664, = 1.25789, 1)4 = -0.51635,176 = 53.8739, c = 0.493339} {<}), 
{vo = 0.998367,172 = 1.35365,174 = -1.83773,176 = 55.4401, c = 0.494999} (+) and 
{i7o = 0.997691, 172 = 1.81431, 174 = 4.97792, i76 = 20.8538, c = 0.493152} (□). 

Figure 7: Umax versus gravity g for four different channel widths L = 32, 48, 64, 128 
(from bottom to top) and for e = 0.02. The lines are least-squares linear fit to the 
numerical data. Four lines cut the vertical axis at non-zero positions which are very 
close to each other. 

Figure 8: (a) The width of higher densities (than average density) relative to the 
channel width A/L are plotted versus the gravity g for different channel widths and 
different dissipation parameters. A: {e = 0.01, L = 32}, B: {e = 0.01, L = 128}, 
C: {e = 0.02, L = 48}, D: {e = 0.02, L = 64}, E: {e = 0.02, L = 128}, F: 
{e = 0.04, L = 48}, G: {e = 0.04, L = 64}. (b) Data collapse of (a): A/L are 
plotted versus g'^L^ /e with yU = 1.4 and v = 0.5. 

Figure 9: a versus L^/^ for e = 0.01 {<}) and e = 0.02 (+). 

Figure 10: Data collapse of Fig. 5: 7 versus F = 10~^g^^'^L^^^'^e. Symbols are the 
same as those in Fig. 5. Solid line is the best fit to the numerical data using Eq. (^). 
The logarithm of horizontal axis is made for clarity of displaying data. The fitting 
parameters in Eq. (|2|) are d = 0.0097748, e = -0.350605 and h = 0.00268941. 
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